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INTRODUCTION 


The objectives of NASA Grant #NAG-1-158, "Response of Composite 
Plates Subjected to Acoustic Loading" were to investigate numerical 
methodology for the determination of narrowband response in the 
geometrically nonlinear regime, to determine response 
characteristics for geometrically nonlinear plates subjected to 
random loading and to compare the predictions with experiments to 
be performed at NASA Langley Research Center. The first two 
objectives were met. The response of composite plates subjected 
to both narrowband and broadband excitation have been studied and 
much useful information has been obtained (as will be discussed in 
this report and is documented in two AIAA papers [1,2]). The third 
objective could not be met due to lack of experimental success and 
the decision to discontinue this research on the part of the 
sponsors. A brief summary of the experimental efforts can be found 
in Appendix A of this report. 

The problem of determining the response of composite panels 
subjected to high intensity acoustic loading is of extreme 
importance for future flight applications involving high speed, 
high altitude flight of structures which contain composite members. 
Several reports and papers have indicated that future flight 
applications are likely to require composite panels to perform 
under acoustic loading conditions which will produce geometrically 


nonlinear responses [3,4,5]. Even though the application is 
crucial, little robust analysis methodology exists which can 
determine the response characteristics in this regime. This report 
addresses a methodology approach which can handle a wide variety 
of applications in a robust, self consistent manner. 

The problem of the geometrically nonlinear response of composite 
panels subjected to narrowband, sinusoidal loading was studied 
initially. It is well known that flat plates exhibit response 
characteristics which are multivalued in the frequency domain. 
This phenomena manifests itself in jump characteristics which can 
be observed experimentally. A numerical approach to the solution 
of this problem was developed and reported previously [1]. This 
work demonstrated that a time domain integration approach with a 
phased damping model and an energy convergence criterion could 
predict the solution quite accurately. In addition, a numerical 
method was proposed and tested for the prediction of jump 
phenomena. This approach did not require the advocation of an ad 
hoc stability criteria. As was expected, the differential equation 
contained all relevant information for the prediction of the jump 
frequencies. It was demonstrated that the nonphysical solution was 
numerically unstable. The typical hard spring response was 
accurately predicted and the metastable transition response 
characteristics were predicted. This paper is contained in 
Appendix B which details this work. 


The problem of the response of a composite plate subjected to 


random loading was studied next. Several different approaches have 
been proposed over the years for the analysis of nonlinear systems 
subjected to random time histories. These include exact (Fokker- 
Planck) equation solutions [6,7], perturbation approaches [8,9], 
stochastic linearization [10,11] and time domain (or Monte Carlo) 
approaches [12,13]. The Fokker-Planck equation solutions and 
perturbation approaches are discussed in the referenced literature. 
Neither methodology, however, can be employed for general loading 
and arbitrary degrees of nonlinearity. These approaches, 
therefore, were not pursued. 

The method of stochastic linearization is based on the assumption 
of a known response density function. While most authors employ 
a Gaussian closure, current work is focusing on non-Gaussian 
closures. The basic approach seeks an approximate solution to the 
nonlinear problem by proposing a "equivalent" linear problem which 
exhibits similar response statistics. While this method is 
analytically attractive since the solution of the equivalent linear 
system is trivial, it requires knowledge of the response density 
function and can only minimize error with respect to a few response 
statistics (most commonly the RMS response). As is shown in [2], 
stochastic linearization produces good results for the response 
statistics for which the error is minimized but does not predict 
the correct spectral characteristics. In addition, the approach 
assumes that the spectral response width is independent of the 
nonlinearity (the peak response frequency, however, is dependent 
on the nonlinearity) which is demonstrated to be incorrect. 



The time domain or Monte Carlo simulation approach appears to be 
the most rigorous and robust approach to the solution of the 
problem of determining the response characteristics of a composite 
plate undergoing geometrically large deflections due to random 
excitation. This approach was compared with the method of 
stochastic linearization and with the exact solution for white 
noise excitation in [2]. The exact response density function 
and response statistics were predicted by the time domain 
simulation method. Much discrepancy exists# however# between the 
exact (and numerical time domain) solution and the stochastic 
linearization solution. The details from [2] are given in Appendix 
C. 


This report will discuss the results obtained from varying the 
degree of geometric nonlinearity# damping and loading amplitude. 
Special focus is placed on the characteristics and coupling of 
these effects. All solutions were generated for white noise 
loading (where the numerically predicted response density matches 
the exact solution) and clamped boundary condition (which cause 
the most severe response and fatigue characteristics). 


TIME DOMAIN SIMULATION APPROACH 


For a composite panel subjected to a uniform pressure loading with 


random time history, the single mode response reduces the governing 
system to the well known Duffing equation given by 
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For comparative purposes, the hard spring coefficient was taken 
for the composite panel discussed in Appendix C. The nonlinearity 
parameter was varied to study the effects of increased and 
decreased nonlinearity. A value of unity corresponds to the 
results in Appendix C. In addition, soft spring characteristics 
could be studied by employing negative nonlinearity parameters. 

The Duffing equation was integrated using the implicit Newmark 
integrator as discussed in Appendix C. During the course of this 
research, several integrators were investigated: explicit second 
order central differences, explicit Runge-Kutta algorithms (both 
fourth and sixth order) and two of the implicit Gear integrators 
(specially designed for stiff systems). These algorithms are 



detailed in numerical analysis text books (e.g. [14]) and will not 
be repeated here. 

Integration with explicit algorithms was unsuccessful. The time 
step size had to be chosen prohibitively small to accommodate 
random loading and maintain stability. Adaptive measures appear 
to be promising, however, due to the nonlinearity, they must be 
developed empirically and would depend on all involved parameters. 
This approach was demonstrated for one response but was still 
computationally expensive compared with the Newmark approach. 
Increasing the order of the integrator did not improve the 
situation as the required time step size was governed by stability 
and not accuracy. The stable time step size appeared to produce 
response statistics and response density functions which were 
within 1% of the exact solution. Explicit integration was not 
pursued further. 

The Gear integrators were tested on the problem also. These 
integrators were developed for integration of stiff equations 
(mainly for application in the field of chemical kinetics). For 
first order stiff systems (which exhibit exponential 
characteristics), these integrators have been very successful. 
For problems with time harmonic characteristics (such as vibration 
problems), they have not been very successful. Both second and 
fourth order Gear integrators were tried for three of the parameter 
sets studied. The convergent time step size was much smaller 
(approximately 20%) than required by the Newmark integrator. In 



addition, the resulting nonlinear algebraic equation was more 
difficult to solve using Newton's method. This was especially true 
for balanced nonlinearity (where the linear and RMS nonlinear 
stiffness were approximately equal). The time step had to be 
reduced almost an order of magnitude before Newton's method would 
converge from time step to time step. The Gear integrators, 
therefore, were abandoned. 

The Newmark integrator proved to be a stable and accurate performer 
for this application. The convergence and staggered damping 
algorithm employed are discussed in Appendix C. The convergence 
characteristics seemed quite insensitive to the degree of 
nonlinearity and damping (for the range studied). The adaptive 
time stepping discussed in Appendix C was also quite insensitive 
to changes in damping and nonlinearity. Since this algorithm was 
developed empirically, a more optimal approach should be pursued 
for application to regimes with significantly different 
characteristics. The parameters chosen produced the exact response 
density (within 1%) for the most severe cases (e.g. highest degree 
of nonlinearity, both extremes in damping and highest RMS 
response) . 

To investigate the effects of varying nonlinearity and damping, 
the time histories were generated by integrating the Duffing 
equation to a steady state and continuing until 1,000,000 samples 
were obtained (as detailed in Appendix C). The RMS response was 
calculated for each case as well as the response spectral density 



(often called the Power Spectral Density or PSD). Since the 
response density was symmetric, only even statistical moments 
contain useful information for white noise loading. For other 
loading distributions, however, odd order statistical moments may 
prove informative. Due to time limitations, only the RMS moment 
was analyzed. This is further discussed in the last section. 

A modified algorithm for Spectral Density estimation was introduced 
to improve the accuracy of the calculation. An overlapping window 
approach was employed which minimized variance significantly. For 
linear problems where the exact PSD is known, the overlapping 
approach reduced variance by almost an order of magnitude. Details 
of this algorithm are given in Appendix C. This approach was 
employed for all results to be presented. 


RESULTS AND DISCUSSION 

For the purposes of this study, three different nonlinearity 
parameters (x s 0.5, 1.0, 1.5) , five different damping factors 
( £ = 0.02 , 0.04 , 0.06 , 0.08 , 0.10 ) and six different 

excitation amplitudes ( 70, 80, 90, 100, 110, 120 dB ) were 
investigated. For each combination of the parameters, the RMS 
response amplitude and Response Spectral Density (RSD) were 
calculated. As previously discussed, the response density function 
is non-Gaussian and known for this special case of white noise (see 
Appendix C) . 



For the nonlinearity parameter of 0.5 (the least nonlinearity 
studied), the RMS response is shown in Figure 1. For all damping 
levels studied (up to 10% of critical), the nonlinearity is evident 
in the RMS response by 110 dB loading level. At the low end of the 
damping, nonlinearity is evident at 100 dB. Though not evident on 
the logarithmic scale, there is measurable (greater than 10%) 
deviation from linearity at 90 dB. 

Figures 2 and 3 show the RMS response for the nonlinearity 
parameter equal to 1.0 and 1.5 . At 1.0, the deviation from the 
linear RMS curve is larger than at 0.5. This is the expected 
result for increasing nonlinearity. It is interesting to note that 
even for loading levels of 90 dB (which is not terribly severe), 
this typical panel (described in Appendix C) demonstrates 
nonlinearity which should not be ignored. As will be evident 
subsequently, this effect can be seen more clearly in the shift of 

the spectral peaks. At a nonlinearity parameter of 1.5, the 

\ 

nonlinearity is evident on the plot at 100 dB. Since the increased 
nonlinearity also increases the stiffness of the plate, the 
amplitude is lower and, therefore, the deviation from the linear 
line is less. This is not to be confused with a lesser degree of 
nonlinearity. 

Figures 4, 5 and 6 show the RMS response at constant damping 
factors of 0.02, 0.06 and 0.10 as a function of the nonlinearity 
parameter. The larger amplitude is predicted at the lower 


nonlinearity parameter since the total stiffness is a monotonically 
increasing function of the nonlinearity parameter. As the damping 
level increases/ the nonlinearity parameter is not as critical a 
factor as the increased damping diminishes the RMS response and 
hence the nonlinear stiffness contribution. As is evident in 
Figure 6, even at a large damping ratio of 10% critical, 
nonlinearity is evident (even on a logarithmic scale) at 100 dB 
which may be evident in many applications. 

Figure 7 shows the RSD as a function of loading amplitude for a 
constant damping ratio of 2% critical and a nonlinearity parameter 
of 1.0 . As the amplitude increases, the spectral peak near the 
linear resonance frequency shifts to a larger frequency which is 
typical of a hard spring nonlinear response. For a loading level 
of 90 dB, the peak has shifted on the order of 10%. It is 
interesting to note that as the peak shifts, it also broadens. 
Spectral broadening has been noticed experimentally by several 
researchers (as discussed in Appendix C). Often, spectral 
broadening is attributed solely to damping phenomena (often modeled 
as nonlinear damping effects). Nonlinear stiffness effects have 
not been considered a mechanism for spectral broadening previously. 

Figures 8 and 9 show the RSD at nonlinearity parameter values of 
0.5 and 1.5 . For the smaller nonlinearity parameter, the RSD is 
larger than for the higher values. This is due to the fact that 
increased nonlinearity parameter increases the total stiffness of 
the oscillator. Increased stiffness decreases the response. 


suppressing the RSD. Each of these plots demonstrates shifting of 
the spectral peaks for loading levels of 90 dB and larger. It is 
interesting to note, that for a loading level of 120 dB, no peaks 
of any form can be discerned. The RSD exhibits a gradual rise from 
the zero frequency value over the range plotted. If the 
calculation were performed to higher frequencies, the RSD is 
expected to gradually decrease past a limiting value of frequency. 
The numerical results generated are valid out to a frequency of 
1000 hz. The RSD calculations are valid out to 500 hz. This 
phenomena was not investigated further in this study. 

Figures 10, 11, 12 and 13 show the RSD at a fixed value of the 
nonlinearity parameter (1.0) and damping levels of 4%, 6%, 8% and 
10% of critical (respectively). As the damping increases, more 
spectral broadening is observed. In addition, the peaks shift less 
at higher damping levels. As damping increases, the spectral 
broadening due to damping is evidenced at low loading levels. At 
high loading levels, however, due to the decreased nonlinear 
stiffness term, the total spectral broadening decreases due to the 
decreased nonlinear stiffness and, therefore, decreased broadening 
due to nonlinear stiffness (this will be more evident when the data 
is shown at constant amplitude and nonlinearity parameter). 

Figures 14 and 15 show the RSD for a damping level of 10% critical 
with nonlinearity parameters of 0.5 and 1.5. For increased 
nonlinearity parameter, the peaks broaden due the increased 
nonlinear stiffness. Even though the response amplitude is 


decreased at the larger nonlinearity parameter (due to the larger 
total stiffness), the spectral broadening is greater. This would 
suggest that the percentage of stiffness due to the nonlinear term 
is a governing parameter for the degree of spectral broadening. 

Figures 16 and 17 show the RSD predicted for a nonlinearity 
parameter of 1.0 and a driving amplitude of 70 and 80 dB. As 
expected, the response is linear. At larger damping ratios the 
peak broadens, however, the frequency of maximum RSD remains the 
same (i.e. no spectral shift is observed). Figures 18 and 19 show 
the RSD at a driving amplitude of 80 dB for nonlinearity parameters 
of 0.5 and 1.5 . As expected, since the response is linear at this 
point, little difference is observed. 

Figures 20, 21 and 22 show the RSD at 90, 100 and 110 dB with the 
nonlinearity parameter equal to 1.0 . As the amplitude increases, 
the peaks shift and broaden at all damping levels. It is clearly 
evident in Figures 21 and 22 that less peak shift and less peak 
shift is evident at higher damping levels. In other words, at a 
constant amplitude and nonlinearity parameter, the spectral peaks 
get narrower at higher damping levels since the additional damping 
decreases the amplitude which decreases the nonlinear stiffness and 
decreases the spectral broadening due to nonlinear stiffness 
effects. 

Figures 23 and 24 show the RSD at 110 dB for nonlinearity 
parameters of 0.5 and 1.5. All exhibit the same decrease in 



spectral width with increasing damping. It is dramatically 
demonstrated in Figure 24. Figure 25 shows the RSD at 120 dB for 
the nonlinearity parameter equal to 1.0 . For large damping, a 
round, broad peak region can be identified. As the damping 
decreases, however, the peak broadens to the point of losing 
distinction. At large amplitude loading, the RSD exhibits wide, 
flat spectral characteristics. 

Figure 26 shows the RSD at 90 dB with 2% critical damping as a 
function of the nonlinearity parameter. As mentioned previously , 
the larger nonlinearity parameter shifts the peak to a greater 
degree than for smaller nonlinearity parameters. This is slightly 
evident in Figure 26. For the same amplitude. Figures 27 and 28 
show the RSD for 6% and 10% critical damping levels. The increased 
damping suppresses the response and lessens the nonlinear stiffness 
effects. At these damping levels, therefore, the effect of the 
nonlinearity parameter is minimal. 

Figure 29 shows the RSD at a loading amplitude of 100 dB and a 
damping level of 6% critical. While not evident at 90 dB, the 
nonlinear stiffness parameter shifts, broadens and flattens the 
peak some at 100 dB for this damping level. Figure 30 shows the 
RSD at 100 dB for 10% critical damping. There is still little 
influence due to the nonlinearity parameter under these conditions. 


Figure 31 shows the RSD at 110 dB with 2% critical damping 



Increasing the nonlinearity parameter significantly shifts, 
broadens and flattens the RSD under these conditions. Figures 32 
and 33 show the RSD at 110 dB for 6% and 10% critical damping. 
While the nonlinearity parameter shifts broadens and flattens the 
response in these figures also, the degree is diminished as a 
function of increased damping. At 110 dB with 10% critical 
damping, the spectral widths, while showing some broadening due to 
nonlinearity parameter, are not very dramatically changed. This 
indicates that there may be a crossover point (for constant 
amplitude loading) where the damping ceases to dominate the 
spectral width and where the nonlinear stiffness effects become 
the dominant broadening mechanism. 


CONCLUDING REMARKS 


The work under NASA Grant #NAG 1-158, "Response of Composite Plates 
Subjected to Acoustic Loading" produced many significant and 
insightful results which should enhance the understanding of 
randomly excited nonlinear oscillators and, ultimately, lead to 
accurate prediction of sonic fatigue characteristics. During the 
two funded years of this research grant (and the six month no cost 
extension), two papers were presented and published with the AIAA. 
In addition, a third paper incorporating all the results presented 
in this report is under preparation for submission to an 
appropriate journal. Upon completion, this paper will be forwarded 
to the NASA Grant Monitors, Dr. John Mixson and Mr. Carl Rucker. 


The results of this research demonstrated that nonlinear stiffness 
effects cause spectral peaks to broaden and flatten as well as to 
shift the frequency of maximum response. This discovery was 
unreported prior to this work. In addition/ it was demonstrated 
that increased damping, while broadening spectral peaks in the 
linear stiffness range, narrows spectral peaks and reduces spectral 
shift where nonlinear stiffness effects are large. In addition, 
as nonlinear effects become more dominant, peak definition is 
tending to total obfuscation. 

The results of this study demonstrates that nonlinear stiffness 
effects may play an important role in the response of composite 
panels even relatively low loading levels. For the panel discussed 
in Appendix C, the nonlinear effects are evident at 90 dB and are 
significant at 100 dB. Many applications currently see acoustic 
loading of these intensities. Consideration of nonlinear effects, 
therefore, may be more important than was previously considered. 

Of important consequence is the fact that nonlinear stiffness 
effects are a significant mechanism of spectral broadening. This 
study demonstrates that not only are the nonlinear stiffness 
contributions significant but that they are coupled to the effects 
of linear damping. While not studied in this work, it is a logical 
conjecture that nonlinear damping mechanisms, if relevant will have 
effects which are coupled to nonlinear stiffness effects. It may 
be unrealistic in application, therefore, to uncouple or ignore 



nonlinear effects even in relatively mild loading regimes. 

For the important application of sonic fatigue life prediction, 
the accurate understanding of panel response is critical, but not 
sufficient knowledge. Just as nonlinear effects play a significant 
role in the response characteristics of the deformation, they may 
have separate and sizeable impact on fatigue life characteristics. 
Most fatigue damage theories in vogue either explicitly or 
implicitly assume the applicability of modal superposition and the 
validity of the Fourier transformation. This has by no means been 
demonstrated for the nonlinear oscillator. The results of this 
study demonstrate, for the white noise loading regime and the 
simple nonlinearity considered, that the resulting deformation 
process are ergodic and stationary. This does not, however, 
indicate that the mapping from time domain to frequency domain (and 
the inverse) are correct or unique. Indeed, the results for 
sinusoidal loading indicate that the mapping is probably not 
unique. 

To fully understand even the single degree of freedom oscillator, 
significant additional research is required and must be performed 
in the time domain. Loading spectra other than white noise must 
be investigated, coupling effects between stiffness and damping 
mechanisms must be studied and cross correlation between field 
quantities (e.g. velocity and deflection) must be delineated. The 
current study, while yielding important insight into the problem, 
has just scratched the surface of understanding the response of 


nonlinear oscillators subjected to random loading 
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APPENDIX A: COMPARISON WITH EXPERIMENTS 



SUMMARY OP EXPERIENCE WITH THE LANGLEY DATA 


One of the original goals of this study was to compare the 
predicted response characteristics with tests to be performed under 
the direction of Mr. Carl Rucker and Dr. John Mixson at the NASA 
Langley research Center. Tests were to be run on both aluminum and 
graphite-epoxy panels. These panels were to be secured in a 
circular aperture (rigid) to insure clamped edge conditions. 

Several tests were performed. The results can be obtained from 
Mr. Rucker. At The George Washington University, we reduced much 
of the data obtained. Preliminary data indicated that white noise 
conditions (i.e. flat loading spectra) and uniform load 
distribution across the panel were not well produced in the lab. 
In addition, the aluminum panel results indicated that the boundary 
conditions were not rigidly clamped edges. It also did not appear 
possible to drive the aluminum panel into the nonlinear region. 
For these reasons, a joint decision was made to not compare 
numerical results with the data for the aluminum panels. It was 
hoped the composite panels would show better boundary response and 
be more easily excited into the nonlinear regime. 

The composite panel results demonstrated additional problems. Due 
to thermal conditions in the laboratory, the panel buckled and 
deformed during the history. This distortion manifested itself in 
soft spring nonlinear behavior for some test runs and hard spring 
on others. It was impossible to determine the appropriate 


parameters for this panel even for a given test run. Modeling of 
initial curvature effects was beyond the scope of this study. In 
addition, it was not clear that the buckling was at all 
geometrically regular or that the boundary conditions were any 
closer to perfectly clamped. A final difficulty was in trying to 
determine an equivalent damping parameter. Repeating tests on the 
panel during a single day produced different results indicating a 
history dependence obviously absent from the analytical model. It 
was jointly determined that comparison with this data would not be 
fruitful. 

In the real world, panels used in aircraft and aerospace structures 
are apt to demonstrate many of the same phenomena exhibited by the 
test panels at NASA Langley. The understanding of these effects, 
however, must be determined after understanding the ideal flat 
plate with known (if not perfectly clamped) and quantifiable 
damping. Response and fatigue models must be compared with a 
reliable data base generated from panels and loading conditions 
which are as close to the theoretical assumptions as possible. For 
understanding of composite panel response to acoustic loading and 
for the prediction of sonic fatigue life, significant testing is 
required under very basic conditions. These tests must be 
performed prior to the introduction of other effects (such as 
thermal expansion, thermal gradients, buckling and curvature 
effects) or it will be impossible to discern the applicability of 
any theoretical or numerical modeling. Correlation with this type 
of complicated data prior to establishing the basic theoretical 


modeling will reduce to an exercise in parameter fitting which can 
not be extrapolated to application. 
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Abstract 

The response of composite panels subjected to random pressure loads 
large enough to cause geometrically nonlinear responses is studied. A 
time domain simulation is employed to solve the equations of motion. 
A staged damping approach is employed to remove transients 

efficiently. An adaptive time stepping algorithm is employedto 
minimize intermittent transients. A modified algorithm for the 
prediction of response spectral density is presented which predicts 
smooth spectral peaks for discrete time histories. 

Results are presented for a number of input pressure levels and 
damping coefficients. Response distributions are calculated and 
compared with the analytical solution of the Fokker-Planck equations. 
RMS response is reported as a function of input pressure level and 
damping coefficient. Spectral densities are calculated for a number of 
the examples. The results show excellent agreement with the 
Fokker-Planck solution. The RMS predictions and the predicted 
distributions are extremely accurate. The spectral density 

calculations demonstrate that the response peaks broaden as and 
flatten as a function of input pressure level (at a constant damping 
level) . These results indicate that the nonlinear stiffness of the 
panel can account for some spectral broadening of the response. 


Introduction 

Over the past twenty five years, many authors have studied the problem 
of obtaining an understanding of how a flexural member responds to 
random, time dependent loading. If the member response is linear, 
exact methodologies exist to address this problem (see, for example, 
[1]). Often, however, the member exhibits nonlinear characteristics 
due either to nonlinear material behavior (e.g. plasticity or 
viscoplasticity) or "large" deformations. In either of these 
situations, the problem becomes much more difficult. 

Of interest in aerospace applications is the problem of determining 
the response of flexural members subjected to random loading in the 
regime where geometric nonlinearities occur. Many current and future 
applications contain critical structural elements which exhibit 
response magnitudes which necessitate nonlinear analysis (see, for 
example, [2,3,4]). A clear understanding of the response of these 
members is critical to the success and safety of current and future 
aerospace vehicle design. 

Several different approaches have been proposed over the years for the 
analysis of nonlinear systems subjected to random time histories. The 
most commonly employed methodologies are the Fokker- Planck equation 
solutions [5,6], perturbation methods [7,8], stochastic linearization 
techniques [9,10] and time domain methodologies (or Monte Carlo 
approaches [11,12]). Each of these methods have their proponents and 
have been employed for a variety of problems. Each, however, have 
some inherent problems or open questions. 


The Fokker- Planck equation solutions are to be employed whenever 



possible as they yield exact solutions for the response probabilities 
and statistics. Unfortunately, the input spectra which can be 
analyzed is extremely limited. The solutions are best employed for 
comparison with approximate methodolgies as a test problem. 

Perturbation approaches are employed for a wide class of nonlinear 
problems in mathematics and mechanics. Since a variety of techniques 
can be devised around the same concept it is difficult to generalize 
as to the range of applicability of the approach. It is fair, 
however, to note that the method is based on the assumption that the 
nonlinearity is small relative to the total characteristic of the 
system. This methodology, therefore, is best suited to weakly 
nonlinear systems. 

Over the years, many methods of stochastic linearization have been 
suggested for the solution of nonlinear oscillator problems. In 
general, peak response frequencies and RMS (root mean square) response 
characteristics are predicted to within engineering satisfaction for 
moderately nonlinear systems. The methodolgy, however, has been 
applied to a wide variety of problems in which the response 
probability density function and the driving probability density 
function are quite different. In these problems, stochastic 
linearization can not yield meaningful spectral information even if 
the peak response frequencies and RMS values are reasonable. 

Time domain simulation or Monte-Carlo approaches have been employed 
for the solution of nonlinear systems for many years. While the 
methodology is straightforward, several problems exist. For systems 
with small structural damping, considerable time may be required to 



damp out the transients in the system. Also if the damping is light, 
errors in the solution can excite pseudo -transients in the response 
which can cause local errors in the response. In addition, depending 
on the time integration approach chosen and the stability limit, large 
numbers of time steps are often required to obtain statistically valid 
response distributions [13,14]. The most critical problem with time 
domain solutions, however, is the difficulty in extracting spectral 
information from discrete time histories. 

In this work, a single oscillator subjected to wide band white noise 
is studied by time domain simulation. The Newmark time integrator is 
employed for reponse prediction [15]. A graded damping approach is 
employed which minimizes the time required to remove the transients in 
the system without affecting the long time response [16] . An adaptive 
time stepping approach is introduced which minimizes the appearance of 
pseudo- transients and improves the local accuracy of the solution. 

A major emphasis in this work is the introduction of an improved 
spectral response calculation approach. The approach employs a 
careful blending of windowing and overlapping data segments which 
reduces (somewhat) computation time while producing a remarkably 
smooth response in the areas of interest. 

The problem addressed is representative of a typical composite panel 
used in application. Response statistics, distribution 
characteristics and spectral densities are predicted. Different 
damping coefficients and driving force levels are studied. Comparison 
is made with other prediction methodologies. 


Problem Description 

Consider a symmetrically laminated composite plate subjected to 
transverse loading. If the plate deflections are assumed to be large 
in the von-Karmen formulation sense, the equations of motion can be 
written in terms of the out of plane displacement, W and a stress 
function F as [17] 


^hW + L^W) - *<F,W) - P(t) - 0 

(1-a) 

L 2 (F) + 1/2 *(W,W) - 0 
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where the differential operators are given by 
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In this formulation, the plate is assumed to be thin and the in-plane 
and rotary inertia contributions have been ignored. The matrix 
components Aij and Dij are dependent only on the material constants . 
If a single space mode for the out of plane displacement, U, is 
assumed to be in the form 

W - - (1 + cos ) (1 + cos -7— ) (3) 

4 a d 

and the stress function, F, is defined accordingly, the equation of 
motion can be reduced to the form 


1> + 2fw J> + wjj* + 01? - 

O O IS 


(4) 


where u> o is the linear angular frequency, 0 is the nonlinear 
stiffness, f is the damping coefficient, m is the total plate mass and 
P(t) is the loading pressure as a function of time. The parameters 
are calculated for a clamped, rectangular panel with uniform pressure 
loading in the transverse direction. The assumed material properties 
are 


*12 
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- 23.7 x 10 6 PSI 

- 1.48 x 10 6 PSI 

- 0.94 x 10 6 PSI 

- 0.30 

1.459 x 10‘ 4 LB - SEC 2 /IN 4 


( 5 ) 


and the plate dimension are 



a - 15 inches 


b - 12 inches (6) 
h - 0.04 inches 

The problem described corresponds to a typical composite panel used in 
application. The damping coefficient is varied in the analysis as 
most application have uncertain damping characteristics. The problem 
described has previously been studied by Mei and Prasad using 
stochastic linearization [18]. Details of the formulation can be 
found in that paper. 

The loading function, P(t) is assumed to be wide band white noise. A 
standard Gaussian white noise generator was implemented [19] to 
produce a flat spectral density. Since a discrete time domain 
solution can only simulate band limited white noise, numerical 
experimentation was employed to determine suitable cutoff frequencies. 
Cutoff frequencies from 10 to 100 times the linear resonance frequency 
were tested. Response statistics were stable for all cutoff 
frequencies studied. Unfortunately, the spectral estimates were quite 
sensitive to cutoff frequency as is well documented. To minimize the 
cutoff frequency effect, all responses employed a cutoff frequency of 
10,000 hz. For practical applications, a lower cutoff frequency could 
probably have been employed, however, spectral calculations may show 
more "noise” as discussed subsequently. 

The Fokker-Planck equations provide an exact solution for this problem 
[6]. For the general oscillator of the form 


X + 0X + F(X) - f(t) 


(7) 



where 


F(X) - «>* [X + eg(X) ] 
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The response distribution function is given as 

2 
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where 


G(X) - 


g(cr)da 


( 10 ) 


For the problem of equation 5, therefore, the distribution is 


P(tf) - C exp ^<1“ + x 2^ 4 ^] 


( 11 ) 


C, X^ , 12"* constant 


This will be compared with the time domain solution to demonstrate the 
accuracy of the response. 

Time. Infiagg&fclgn Approach 

To integrate the equation of motion (equation 5), the Newmark time 
integrator [15] was employed. For the problem in this paper, the 
algorithm takes the form 
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where At is the time step size. This integrator is unconditionally 
stable for the' linear oscillator ( j3 - 0 ). For the nonlinear 
oscillator, however, no stability analysis has been performed to date. 
Experience has shown, however, that accuracy is the limiting concern 
for the Newmark integrator for practical problems [16]. 


The input load is calculated at a discrete number of points determined 
by the cutoff frequency. The data point spacing is given by 
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To solve the equation of motion, therefore, this is the largest time 
step permissable . In practice, choosing the time of equation 11 will 
yield a solution with intermittent transients. Numerical 
experimentation demonstrated that the time step required to solve the 
equation of motion without encountering intermittent transients was on 
the order of 

At - (14) 

for the most severe case studied. To produce a statistically 
significant number of points, however, would have required a large 
number of time steps. To reduce this requirement, an adaptive 
technique was introduced. Assume the solution is known at a time T 
with corresponding driving function Q(T) . It is desired to proceed 



from T to T + Ar in H equal time steps. The driving function at T 
+ Ar is given as Q(T+ Ar). If Qmax and Qmin are the maximum and 
minimum values the driving function exhibit, the number of time steps 
is chosen as 


M - 2 + 10 * 


KT + Ar) - Q(T) 
^max Snin 


(15) 


This worked well for all problems studied. The choice of the 
multiplicative factor in equation 13 should, however, be a function of 
the damping coefficient and the driving function level instead of 
being held constant. Since this was chosen for the most severe case, 
the results should be good for all cases. More work could lead to an 
optimization of this approach, however. 


In the time domain solution approach, damping is required to produce a 
steady state solution. While it can be argued that damping will be 
exhibited by any real structure, that damping is often light. To 
remove the transients from the system may require significant amounts 
of computation. To accelerate removal of the transients, a staged 
approach was employed. To start the solution, initial response and 
velocity was assumed to be zero. The damping coefficient was 
multiplied by 1000 and the solution was generated for 100 * Ar The 
damping was then abruptly reduced to 100 times the true damping level 
and the solution was continued for 100 * Ar . The damping was 
again dropped to 10 times the true damping and the solution continued 
for an additional 100 * Ar. The damping was then set to the true 
(physical level) . The solution was then continued for the desired 
length of time. RMS statistics were accumulated dynamically. When 
both response and velocity values remained constant for time segments 


of 100 times the linear period, the solution was considered to be in 
steady state. The response was then generated to the desired stopping 
point. This was achieved in under 50,000 * A r for all problems 
studied. 


After steady state was achieved, the solution was generated for 
1,000,000 * A r . Samples were taken in Ar increments. These 
points were then used to calculate response statistics and 
distributions. Part of this time trace was also used to calculate 
response spectral density. The approach described is by no means 
optimized. Currently (as described), the parameters are all based on 
numerical experiment. Those employed in this study, however, produce 
very accurate responses with no detectable intermittent transients. 
The cumulative response statistics vary by less than 0.1% over the 
total time history. This methodology could be studied parametrically 
to determine general forms for the parameters involved. 


Spectral PtnsltY Eg£iBtt£l9,n 

A response measure which often yields important information in the 
study of random vibrations is the Spectral Density (often called the 
Power Spectral Density, PSD). The PSD is defined as the normalized 
Fourier Transform of the autocorrelation function. The 
autocorrelation function is defined by 
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and is related to the PSD by 
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Using the theory of Fast Fourier Transforms (FFTs) , the PSD can be 
estimated N/2+1 discrete frequencies equally divided between 0 and the 
cutoff (or Nydquist) frequency by the formulae [20], 
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In practice, this approach leads to much "leakage” which can be solved 
by windowing the data or using a weighted approach to minimize the 
leakage. The algorithm is then modified to 
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The window function , w^ Is chosen to minimize leakage and to have the 
properties of being 0 at the ends of the data segment, 1 at the center 
and be mono tonic between the ends and the center. Generally, windows 
are positive definite functions. The choice of a window function, 
however, Is quite arbitrary. In this work, the Parzen window has been 
chosen as It Is a consistently good performing function for discrete 
oscillator systems [21]. The Parzen window Is given by the formula 
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For N discrete data points separated by the time Increment appropriate 
for the chosen Nydquist frequency, the algorithm presented will yield 
an estimate of the PSD at N/2+1 discrete frequencies. This approach, 
however, yields estimates with horrendous variance for most realistic 
applications. 


To minimize the variance at the discrete estimate frequencies, an 
alternative scheme can be employed. The available data can be broken 
into segments. Each segment can be employed to calculate a PSD 
estimate at M discrete frequencies where 2*M is the segment length. 


The estimates can then be averaged to produce a smoother PSD estimate 
than if data had been used as previously described. If the total 
number of segments sampled is K, then, a total of 2KM points are used. 
It can be shown that this reduces the variance by a factor of K [22] . 
An improvement to this approach involves using overlapping data 
segments. If the total available data is partitioned in samples of 
length M and if K segments are considered, each of length 2M, the 
segment data is chosen such that the segments overlap. In the 
standard approach, the first and second sample make up the first 
segment, the second and third sample make up the second segment, etc. 
This reduces the variance by a factor of approximately 9K/11 for a 
total sample of (K+1)M points [22] . 

For the problem studied in this work, the PSD results showed large 
variance (on the order of 100%) for the linear oscillator using the 
overlapping segment approach (other methods were worse) . This is 
consistent with other investigator's observations [ 23 ] . To improve 
the PSD estimation, a different overlapping approach was employed. The 
total sample was divided into samples of length J. The first segment 
employed 4 samples: J— 1, J-2, J—3 and J-4. The second segment 
employed samples: J-2, J-3, J— 4, and J-5, etc. For a total of K 
segments, (K+3)*J points were employed. Typical runs for the linear 
oscillator using a Nydquist frequency of 10 times the resonance 
frequency and on the order of 200,000 points, the observed variance 
was reduced to less than 20% in the region near the peak. Specific 
results will be discussed in the next section. 

The approach described above was implemented on an Alliant FX8/2 
computer. Typical runtimes were on the order of one CPU minute for 



200,000 data points. Many numerical experiments were performed both 
on the linear oscillator system and band limited white noise. In 
general, the algorithm performed very consistently over a wide range 
of input signals. Maximum variances were kept below 20% in regions of 
interest and spurious peaks were not observed. White noise results 
were excellent with typical tails on the order of 5% of the Nydquist 
frequency. For a fixed number of data points, the algorithm 
employing three overlapping samples produced significantly better 
estimates for the mean squared response compared with the single 
overlap algorithm (calculated as the discrete integral of the PSD 
using Simpson's 3/8 rule). For single degree of freedom oscillators, 
a prescribed variance tolerance could be met with significantly fewer 
data points using three overlapping segments making this approach 
computationally more attractive. In practice, small tolerances could 
not be met using the single overlap algorithm. 


RgauUa 

The plate described in equations (4,5,6) was studied for many input 
amplitudes and damping coefficients. Figure 1 is a plot of the RMS 
response amplitude 



N 


j-1 (N^I) 


(23) 


as a function of the driving Sound Spectrum Level (SSL) for damping 
coefficients corresponding to 1%, 2%, 4%, 6%, 8% and 10% critical 
(i.e. C - 0.01 , 0.02 , 0.04 , 0.06 , 0.08 , 0.10). The straight 
line represents the results from the linear theory with f-0.01 . As 
is well known, the response predicted by the linear theory is greater 
than that predicted by the nonlinear theory for a given damping 
coefficient. In addition, the increase in damping causes a decrease 



in the RMS response as is expected. The prediction are in good 
agreement with the predictions using the Method of Equivalent 
Linearization (MEL) reported by Mei and Prasad [18 ]. Maximum 
difference between the two methods was 4%. 

To further investigate the accuracy of the results obtained by time 
domain simulation, the probability density function was calculated 
from the response sample. The range of responses ( -3 < V < 3 ) was 
divided into 200 bins equally distributed. The histogram was then 
integrated and normalized numerically so the total probability 
integral equaled 1. The points at the center of each bin where then 
plotted and connected. A total of 1,000,000 sample points were 
employed. The results are shown in Figure 2 for an input SSL of 110 
DB and a damping coefficient of f -0.01. In addition, the exact 
distribution of equation 9 is plotted along with the predicted 
Gaussian response of the MEL. The results from the time domain 
simulation are in remarkable agreement with the exact solution as is 
shown. The Gaussian prediction, however, is quite different. The MEL 
predictions heavily weight the response toward the zero mean and 
overestimate the tail distribution. The time domain results, in 
contrast, accurately predict the broader and flatter distribution of 
the exact solution. The exact solution is modeled quite well over the 
whole distribution range, including the tails. 

Of major interest was the prediction of the response PSD using the 
algorithm described previously. Figure 3 shows the predicted linear 
oscillator PSD for a panel with a damping coefficient of 1% critical 
and a driving input SSL of 70 DB. The response should be a linear 
damped oscillator respone with exact solution in the form 



(24) 
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where is the input spectral density corresponding to 70 DB SSL. 
The predicted PSD is in good agreement with the exact solution for the 
range plotted. The PSD varies by almost three orders of magnitude in 
the range studied near the peak. The response PSD peak and shape are 
well predicted. As will be discussed subsequently, PSD tails tend to 
be over predicted and quite noisy. For physical interest, however, 
this is not of major consequence since this happens when PSD values 
are several orders of magnitude below peak values. 

Of interest is the effect of increasing driving SSL on the response 
PSD. The MEL results of Mei and Prasad indicated that increasing SSL 
did not effect the spectral characteristics of the oscillator. MEL 
predicts a shift in the resonance frequency and in the magnitude of 
the PSD, but not in the width or shape of the distribution. Figure 4 
shows the PSD from the time domain simulation and from the MEL for a 
driving SSL of 110 DB with 1% critical damping. The two methods 
predict quite different response PSDs . The time domain results 
demonstrate significant broadening and flattening of the spectral 
density. The characteristic linear oscillator sharp peak is no longer 
evident. The MEL results predict the expected linear sharp peak. The 
MEL results predict a much larger peak value also. It is important to 
notice that the PSD predictions from the time domain simulation show 
relatively small variances and consistent response patterns. The 
difference between the two methods does not appear to be due to the 


numerical approach employed. 


Figure 5 is a plot of the response PSD for the nonlinear oscillator 
with 1% critical damping at 80, 90, 100 and 110 DB SSL driving levels. 

The 80 DB results are fairly linear in form with a resonance 
frequency of around 94 hz. (the linear resonance frequency is 92.9 
hz.). At 90 DB, the PSD results have shifted slightly (peak response 
is now around 100 hz.) and the peak has broadened slightly. At 100 
DB, a larger shift is observed and the response has broadened and 
flattened considerably. A discernable resonance peak is no longer 
observable. As previously discussed, at 110 DB significant broadening 
and flattening has taken place. All results in this figure are for a 
constant damping ratio (f-0.01). The PSD results generally are 
considered by the author to be accurate. It is important to mention, 
however, that the tails may be overestimated. Studies with linear 
oscillators (with known exact solutions) demonstrated overestimation 
of tail predictions but only when the response was several orders of 
magnitude less than the peak. 

Also of interest is the effect of damping on the predicted response 
PSDs. Figure 6 is a plot of the predicted response PSDs at 1%, 2%, 
4%, 6% and 10% critical damping coefficients. As Is expected from the 
linear theory, the peaks diminish and flatten with increasing damping 
(linear, viscous, damping only is considered in this study). Damping 
only effects the peak responses and the tails are insensitive to 
damping as expected from the linear theory. It is evident from Figure 
1 that at 80 DB SSL, the system should behave linearly. This is 


confirmed. 



Figure 7 is a plot of the predicted response PSDs at 1%, 2%, 4%, 6% 
and 10% critical damping coefficients. In contrast to the linear 
oscillator, the increase in damping decreases the spectral width and 
sharpens the response PSD. The high frequency side of the PSD is 
effected much more than the lower frequency side. The frequency at 
maximum PSD moves back toward the linear resonance frequency, in 
effect, linearizing the system. This is not totally unexpected. At 
higher damping ratios, the RMS response is lowered, lowering the 
effects on nonlinearity. At 10% critical damping and 110 DB SSL, the 
RMS amplitude is about equal to the RMS amplitude for 1% critical 
damping at 100 DB. The predicted peak frequencies are about the same 
and the predicted maximum PSD values are similar. 

Discussion sod Conclusions 

The results presented demonstrate that the time domain simulation (or 
Monte Carlo) approach provides an accurate numerical methodology for 
the study on nonlinear oscillators. Not only are qualitative 
predictions possible but the exact probability density and response 
statistics are reproducible. The method, therefore, should provide a 
quantitatively accurate tool for studying nonlinear phenomena in other 
loading regimes (i.e. for more complicated loading spectra than white 
noise) . 

Methodology for the prediction of response PSDs was briefly reviewed. 
A new (to the author's knowledge) overlapping scheme is employed which 
reduces prediction variance significantly (but in a quantitatively 
unknown amount) . The results are in good agreement with the exact 
results for linear oscillators. Quantitatively accurate estimates 
have been obtained for several problems with known spectral densities. 



With Che exception of tail regions or within 5- 10% of the Nydquist 
frequency, the results are believed to be accurate for a wide variety 
of problems. 

The spectral density results as a function of increasing SSL were 
studied. At higher SSL, the nonlinearity manifests itself by 
broadening and flattening the PSDs. The resonance frequency (or peak 
frequency) is shifted to higher levels and the broadening is skewed to 
the high frequency side. The drop off is faster for the high 
frequency tail and the response is not symmetric about the resonance 
frequency. It is interesting to note that the BUS response for a 
sinusoidally driven oscillator qualitatively demonstrates the same 
behavior due to hard spring nonlinearities. The resonance peak is 
shifted to higher frequency, the peak is broadened and bent to the 
high frequency side. The broadening and flattening of spectral peaks 
due to hard spring nonlinearity, therefore, may not be totally 
surprising. 

The response results have been compared with those predicted by the 
MEL. MEL predictions agree in RMS responses. PSD predictions, 
however, are quite different. Resonance frequencies (or frequencies 
at peak PSD) are roughly the same, however, at 110 DB, the PSD 
predictions show such a broad resonance that the concept may not be 
applicable. It is not surprising that the RMS response predictions 
are in good agreement as the MEL is based on minimizing the error of 



the RMS predictions. The inherent assumption of a Gaussian response, 
however, is clearly incorrect for large SSL as indicated in Figure 2 . 
This may explain the large difference in PSD predictions between the 
two methods. Since the time domain simulation agrees with the exact 
probability density extremely well, the PSD predictions from the time 
domain are expected to be more reliable than the MEL predictions. 

Increased viscous damping has a linearizing effect on a nonlinear 
oscillator. Since the RMS response is decreased, the nonlinear term 
decreases in importance causing the system to be more linearly. For 
the 110 DB SSL results (shown in Figure 7), the PSD peak is narrowed 
and shifted back toward the linear resonance frequency at higher 
damping ratios. For a linear oscillator, the PSD peak is broadened 
for increasing damping. This broadening, however, is small relative 
to the broadening due to the nonlinear stiffness effects. The 
broadening due to damping, therefore, will probably only be evident 
when the nonlinearity is small or negligeable. 

The results of this study indicate that considerable spectral 
broadening can be attributed to nonlinear stiffness effects alone. 
This may be of significant importance in the prediction of sonic 
fatigue failures. Spectral broadening (and increased probability of 
higher amplitude oscillations) is indicative of a wider peak 
distribution. Sonic fatigue predictions are based on the concept of 
acumulated damage which is a function of the weighted average of the 
instantaneous (or discrete peak) responses. The effect of hard spring 
nonlinearity, therefore, may be more critical than indicated by a 
particular RMS response level. 
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